% For "The Climate Risk Premium" (Lemoine, JAERE, 2020).

% Rather than using analytic expressions, this script adds a unit of time 0 CO2 and
% measures total social cost of carbon as difference in utility.
% Also tests assumption of exogenous CO2 trajectory.



%% Initialize

counter_fine1 = 0;
counter_fine2 = 0;

if ~strcmpi(option_otherloop,'montecarlo_cs') && ~strcmpi(option_otherloop,'montecarlo_damage') % standard case
    
    Paths_perturb.M1 = zeros(Params.montecarlo_draws,max_time,size(nodes_climsens,1),size(nodes_damage,1));
    Paths_perturb.M2 = zeros(Params.montecarlo_draws,max_time,size(nodes_climsens,1),size(nodes_damage,1));
    Paths_perturb.temp = zeros(Params.montecarlo_draws,max_time,size(nodes_climsens,1),size(nodes_damage,1));
    Paths_perturb.cons = zeros(Params.montecarlo_draws,max_time,size(nodes_climsens,1),size(nodes_damage,1));
    Paths_perturb.cons_pc = zeros(Params.montecarlo_draws,max_time,size(nodes_climsens,1),size(nodes_damage,1));
    
    Paths_perturb.M1_exogCO2 = zeros(Params.montecarlo_draws,max_time,size(nodes_climsens,1),size(nodes_damage,1));
    Paths_perturb.M2_exogCO2 = zeros(Params.montecarlo_draws,max_time,size(nodes_climsens,1),size(nodes_damage,1));
    Paths_perturb.temp_exogCO2 = zeros(Params.montecarlo_draws,max_time,size(nodes_climsens,1),size(nodes_damage,1));
    Paths_perturb.cons_exogCO2 = zeros(Params.montecarlo_draws,max_time,size(nodes_climsens,1),size(nodes_damage,1));
    Paths_perturb.cons_pc_exogCO2 = zeros(Params.montecarlo_draws,max_time,size(nodes_climsens,1),size(nodes_damage,1));
            
    if ~strcmpi(option_damagetype,'base')
        Paths_perturb.cons_net = zeros(Params.montecarlo_draws,max_time,size(nodes_climsens,1),size(nodes_damage,1));
        Paths_perturb.cons_net_exogCO2 = zeros(Params.montecarlo_draws,max_time,size(nodes_climsens,1),size(nodes_damage,1));
    end
    
    if option_ocean
        Paths_perturb.oceantemp = zeros(Params.montecarlo_draws,max_time,size(nodes_climsens,1),size(nodes_damage,1));
        Paths_perturb.oceantemp_exogCO2 = zeros(Params.montecarlo_draws,max_time,size(nodes_climsens,1),size(nodes_damage,1));
    end
    
else % are doing Monte Carlo draws from cs feedback or damage distributions for illustrative purposes

    Paths_perturb.M1 = zeros(Params.montecarlo_draws,max_time);
    Paths_perturb.M2 = zeros(Params.montecarlo_draws,max_time);
    Paths_perturb.temp = zeros(Params.montecarlo_draws,max_time);
    Paths_perturb.cons = zeros(Params.montecarlo_draws,max_time);
    Paths_perturb.cons_pc = zeros(Params.montecarlo_draws,max_time);
    
    Paths_perturb.M1_exogCO2 = zeros(Params.montecarlo_draws,max_time);
    Paths_perturb.M2_exogCO2 = zeros(Params.montecarlo_draws,max_time);
    Paths_perturb.temp_exogCO2 = zeros(Params.montecarlo_draws,max_time);
    Paths_perturb.cons_exogCO2 = zeros(Params.montecarlo_draws,max_time);
    Paths_perturb.cons_pc_exogCO2 = zeros(Params.montecarlo_draws,max_time);
    
    if ~strcmpi(option_damagetype,'base')
        Paths_perturb.cons_net = zeros(Params.montecarlo_draws,max_time);
        Paths_perturb.cons_net_exogCO2 = zeros(Params.montecarlo_draws,max_time);
    end
    
    if option_ocean
        Paths_perturb.oceantemp = zeros(Params.montecarlo_draws,max_time);
        Paths_perturb.oceantemp_exogCO2 = zeros(Params.montecarlo_draws,max_time);
    end
    
end


%% Calculate time paths for perturbed path

for index_time = 1:max_time
    
    % Starting time
    if index_time==1
        
        if ~strcmpi(option_otherloop,'montecarlo_cs') && ~strcmpi(option_otherloop,'montecarlo_damage') % standard case
            
            % adjusted initial CO2 value for perturbation to CO2
            Paths_perturb.M1(:,index_time,:,:) = ones(Params.montecarlo_draws,1,size(nodes_climsens,1),size(nodes_damage,1))*(Params.M1_start + Params.perturbation*Params.M1_frac);
            Paths_perturb.M2(:,index_time,:,:) = ones(Params.montecarlo_draws,1,size(nodes_climsens,1),size(nodes_damage,1))*(Params.M2_start + Params.perturbation*Params.M2_frac);
            Paths_perturb.temp(:,index_time,:,:) = ones(Params.montecarlo_draws,1,size(nodes_climsens,1),size(nodes_damage,1))*Params.temp_start;
            Paths_perturb.cons(:,index_time,:,:) = ones(Params.montecarlo_draws,1,size(nodes_climsens,1),size(nodes_damage,1))*Params.cons_start;
            Paths_perturb.cons_pc(:,index_time,:,:) = ones(Params.montecarlo_draws,1,size(nodes_climsens,1),size(nodes_damage,1))*Params.cons_start/Paths.pop(index_time);
            
            Paths_perturb.M1_exogCO2(:,index_time,:,:) = ones(Params.montecarlo_draws,1,size(nodes_climsens,1),size(nodes_damage,1))*(Params.M1_start + Params.perturbation*Params.M1_frac);
            Paths_perturb.M2_exogCO2(:,index_time,:,:) = ones(Params.montecarlo_draws,1,size(nodes_climsens,1),size(nodes_damage,1))*(Params.M2_start + Params.perturbation*Params.M2_frac);
            Paths_perturb.temp_exogCO2(:,index_time,:,:) = ones(Params.montecarlo_draws,1,size(nodes_climsens,1),size(nodes_damage,1))*Params.temp_start;
            Paths_perturb.cons_exogCO2(:,index_time,:,:) = ones(Params.montecarlo_draws,1,size(nodes_climsens,1),size(nodes_damage,1))*Params.cons_start;
            Paths_perturb.cons_pc_exogCO2(:,index_time,:,:) = ones(Params.montecarlo_draws,1,size(nodes_climsens,1),size(nodes_damage,1))*Params.cons_start/Paths.pop(index_time);
            
            if ~strcmpi(option_damagetype,'base')
                Paths_perturb.cons_net(:,index_time,:,:) = bsxfun(@times,ones(Params.montecarlo_draws,1,size(nodes_climsens,1),size(nodes_damage,1)),Params.cons_start*(1-bsxfun(@min,Params.max_pctloss/100,bsxfun(@times,Params.cons_damage,Paths_perturb.temp(:,index_time,:,:).^2))));
                Paths_perturb.cons_pc(:,index_time,:,:) = Paths_perturb.cons_net(:,index_time,:,:)/Paths.pop(index_time);
                Paths_perturb.cons_net_exogCO2(:,index_time,:,:) = bsxfun(@times,ones(Params.montecarlo_draws,1,size(nodes_climsens,1),size(nodes_damage,1)),Params.cons_start*(1-bsxfun(@min,Params.max_pctloss/100,bsxfun(@times,Params.cons_damage,Paths_perturb.temp_exogCO2(:,index_time,:,:).^2))));
                Paths_perturb.cons_pc_exogCO2(:,index_time,:,:) = Paths_perturb.cons_net_exogCO2(:,index_time,:,:)/Paths.pop(index_time);
            end
            
            if option_ocean
                Paths_perturb.oceantemp(:,index_time,:,:) = ones(Params.montecarlo_draws,1,size(nodes_climsens,1),size(nodes_damage,1))*Params.oceantemp_start;
                Paths_perturb.oceantemp_exogCO2(:,index_time,:,:) = ones(Params.montecarlo_draws,1,size(nodes_climsens,1),size(nodes_damage,1))*Params.oceantemp_start;
            end
   
            
        else % are doing Monte Carlo draws from cs or damage distributions for illustrative purposes
            
            Paths_perturb.M1(:,index_time) = ones(Params.montecarlo_draws,1)*(Params.M1_start + Params.perturbation*Params.M1_frac);
            Paths_perturb.M2(:,index_time) = ones(Params.montecarlo_draws,1)*(Params.M2_start + Params.perturbation*Params.M2_frac);
            Paths_perturb.temp(:,index_time) = ones(Params.montecarlo_draws,1)*Params.temp_start;
            Paths_perturb.cons(:,index_time) = ones(Params.montecarlo_draws,1)*Params.cons_start;
            Paths_perturb.cons_pc(:,index_time) = ones(Params.montecarlo_draws,1)*Params.cons_start/Paths.pop(index_time);
            
            Paths_perturb.M1_exogCO2(:,index_time) = ones(Params.montecarlo_draws,1)*(Params.M1_start + Params.perturbation*Params.M1_frac);
            Paths_perturb.M2_exogCO2(:,index_time) = ones(Params.montecarlo_draws,1)*(Params.M2_start + Params.perturbation*Params.M2_frac);
            Paths_perturb.temp_exogCO2(:,index_time) = ones(Params.montecarlo_draws,1)*Params.temp_start;
            Paths_perturb.cons_exogCO2(:,index_time) = ones(Params.montecarlo_draws,1)*Params.cons_start;
            Paths_perturb.cons_pc_exogCO2(:,index_time) = ones(Params.montecarlo_draws,1)*Params.cons_start/Paths.pop(index_time);
            
            if ~strcmpi(option_damagetype,'base')
                Paths_perturb.cons_net(:,index_time) = bsxfun(@times,ones(Params.montecarlo_draws,1),Params.cons_start*(1-bsxfun(@min,Params.max_pctloss/100,bsxfun(@times,Params.cons_damage,Paths_perturb.temp(:,index_time).^2))));
                Paths_perturb.cons_pc(:,index_time) = Paths_perturb.cons_net(:,index_time)/Paths.pop(index_time);
                Paths_perturb.cons_net_exogCO2(:,index_time) = bsxfun(@times,ones(Params.montecarlo_draws,1),Params.cons_start*(1-bsxfun(@min,Params.max_pctloss/100,bsxfun(@times,Params.cons_damage,Paths_perturb.temp_exogCO2(:,index_time).^2))));
                Paths_perturb.cons_pc_exogCO2(:,index_time) = Paths_perturb.cons_net_exogCO2(:,index_time)/Paths.pop(index_time);
            end
            
            if option_ocean
                Paths_perturb.oceantemp(:,index_time) = ones(Params.montecarlo_draws,1)*Params.oceantemp_start;
                Paths_perturb.oceantemp_exogCO2(:,index_time) = ones(Params.montecarlo_draws,1)*Params.oceantemp_start;
            end
            
        end
    end        
    
    % simulate stochastic processes, first for endogenous CO2 case    
    if index_time < max_time
        
        M1_value = Paths_perturb.M1(:,index_time,:,:);
        M2_value = Paths_perturb.M2(:,index_time,:,:);
        temp_value = Paths_perturb.temp(:,index_time,:,:);
        cons_value = Paths_perturb.cons(:,index_time,:,:);
        if option_ocean
            oceantemp_value = Paths_perturb.oceantemp(:,index_time,:,:);
        end
        
        for index_interval = 1:1/Params.timestep
            
            counter_fine1 = counter_fine1 + 1;
            counter_fine = counter_fine1;
            
            % Calculate changes in variables
            run sub_transitions;
            
            % Update values of variables
            M1_value = M1_value + M1change;
            M2_value = M2_value + M2change;
            temp_value = temp_value + tempchange;
            cons_value = cons_value + conschange;
            if option_ocean
                oceantemp_value = oceantemp_value + oceantempchange;
            end
            
            
        end % end stochastic process simulation
        
        cons_value(cons_value>1e300) = 1e300;
        
        % store next period's value
        Paths_perturb.M1(:,index_time+1,:,:) = M1_value;
        Paths_perturb.M2(:,index_time+1,:,:) = M2_value;
        Paths_perturb.temp(:,index_time+1,:,:) = temp_value;
        if option_ocean
            Paths_perturb.oceantemp(:,index_time+1,:,:) = oceantemp_value;
        end
        switch option_damagetype
            case 'base'
                Paths_perturb.cons(:,index_time+1,:,:) = cons_value;
                Paths_perturb.cons_pc(:,index_time+1,:,:) = cons_value/Paths.pop(index_time+1);
            otherwise
                Paths_perturb.cons(:,index_time+1,:,:) = cons_value;
                Paths_perturb.cons_net(:,index_time+1,:,:) = cons_value.*(1-bsxfun(@min,Params.max_pctloss/100,bsxfun(@times,Params.cons_damage,temp_value.^2)));
                Paths_perturb.cons_pc(:,index_time+1,:,:) = Paths_perturb.cons_net(:,index_time+1,:,:)./Paths.pop(index_time+1);
        end
    
    end
      
        
    % simulate stochastic processes, now for exogenous CO2 case
    if index_time < max_time
        
        M1_value = Paths.M1(:,index_time,:,:) + Params.perturbation*Params.M1_frac;  % use original CO2 trajectory, with decayed emission added to it
        M2_value = Paths.M2(:,index_time,:,:) + Params.perturbation*Params.M2_frac*exp(-Params.co2_decay*(index_time-1));  % use original CO2 trajectory, with decayed emission added to it
        temp_value = Paths_perturb.temp_exogCO2(:,index_time,:,:);
        cons_value = Paths_perturb.cons_exogCO2(:,index_time,:,:);
        if option_ocean
            oceantemp_value = Paths_perturb.oceantemp_exogCO2(:,index_time,:,:);
        end
        
        for index_interval = 1:1/Params.timestep
            
            counter_fine2 = counter_fine2 + 1;
            counter_fine = counter_fine2;
            
            % Calculate changes in variables
            run sub_transitions;
            
            % Update values of variables
            M1_value = M1_value + M1change;
            M2_value = M2_value + M2change;
            temp_value = temp_value + tempchange;
            cons_value = cons_value + conschange;
            if option_ocean
                oceantemp_value = oceantemp_value + oceantempchange;
            end
            
            
        end % end stochastic process simulation
        
        cons_value(cons_value>1e300) = 1e300;
        
        % store next period's value
        Paths_perturb.M1_exogCO2(:,index_time+1,:,:) = Paths.M1(:,index_time+1,:,:) + Params.perturbation*Params.M1_frac; % use original CO2 trajectory, with decayed emission added to it
        Paths_perturb.M2_exogCO2(:,index_time+1,:,:) = Paths.M2(:,index_time+1,:,:) + Params.perturbation*Params.M2_frac*exp(-Params.co2_decay*(index_time+1-1));  % use original CO2 trajectory, with decayed emission added to it
        Paths_perturb.temp_exogCO2(:,index_time+1,:,:) = temp_value;
        if option_ocean
            Paths_perturb.oceantemp_exogCO2(:,index_time+1,:,:) = oceantemp_value;
        end
        switch option_damagetype
            case 'base'
                Paths_perturb.cons_exogCO2(:,index_time+1,:,:) = cons_value;
                Paths_perturb.cons_pc_exogCO2(:,index_time+1,:,:) = cons_value/Paths.pop(index_time+1);
            otherwise
                Paths_perturb.cons_exogCO2(:,index_time+1,:,:) = cons_value;
                Paths_perturb.cons_net_exogCO2(:,index_time+1,:,:) = cons_value.*(1-bsxfun(@min,Params.max_pctloss/100,bsxfun(@times,Params.cons_damage,temp_value.^2)));
                Paths_perturb.cons_pc_exogCO2(:,index_time+1,:,:) = Paths_perturb.cons_net(:,index_time+1,:,:)./Paths.pop(index_time+1);
        end
                
        clear co2_value temp_value cons_value;
        clear co2change forcing tempchange conschange;
        
    end
    
    
end % end time looping


%% Report time spent

elapsedtime2 = toc;
disp(['Done looping through time for perturbed path.  Time spent: ' num2str(elapsedtime2/60) ' minutes.']);


%% Utility and scc calculations

%tic; 

disp(['Beginning social cost of carbon calculations.'])

switch option_damagetype
    case 'base'
        Paths.margutility = bsxfun(@times, Paths.cons.^(-Params.power_cons), (Paths.cons - Paths_perturb.cons)/Params.perturbation  ); % simulation x time x climsens_node x damage_node
        Paths.margutility_exogCO2 = bsxfun(@times, Paths.cons.^(-Params.power_cons),  (Paths.cons - Paths_perturb.cons_exogCO2)/Params.perturbation  ); % simulation x time x climsens_node x damage_node
        Paths.consdiff = Paths.cons - Paths_perturb.cons; % simulation x time x climsens_node x damage_node
        Paths.consdiff_exogCO2 = Paths.cons - Paths_perturb.cons_exogCO2; % simulation x time x climsens_node x damage_node
    otherwise        
        Paths.margutility = bsxfun(@times, Paths.cons_net.^(-Params.power_cons), (Paths.cons_net - Paths_perturb.cons_net)/Params.perturbation  ); % simulation x time x climsens_node x damage_node
        Paths.margutility_exogCO2 = bsxfun(@times, Paths.cons_net.^(-Params.power_cons),  (Paths.cons_net - Paths_perturb.cons_net_exogCO2)/Params.perturbation  ); % simulation x time x climsens_node x damage_node
        Paths.consdiff = Paths.cons_net - Paths_perturb.cons_net; % simulation x time x climsens_node x damage_node
        Paths.consdiff_exogCO2 = Paths.cons_net - Paths_perturb.cons_net_exogCO2; % simulation x time x climsens_node x damage_node

end


% To get scc, sum probability-weighted marginal utility change along damage dimension, 
% then probability-weighted sum along climate sensitivity dimension, 
% then average along Monte Carlo simulations, 
% then finally integrate over time

% First taking expectations over damages
margutility_dam = sum( bsxfun(@times,reshape(weights_damage,[1 1 1 size(weights_damage,1)]),Paths.margutility), 4); % simulation x time x climsens_node x 1
margutility_dam_exogCO2 = sum( bsxfun(@times,reshape(weights_damage,[1 1 1 size(weights_damage,1)]),Paths.margutility_exogCO2), 4); % simulation x time x climsens_node x 1
consdiff_dam = sum( bsxfun(@times,reshape(weights_damage,[1 1 1 size(weights_damage,1)]),Paths.consdiff), 4); % simulation x time x climsens_node x 1
consdiff_dam_exogCO2 = sum( bsxfun(@times,reshape(weights_damage,[1 1 1 size(weights_damage,1)]),Paths.consdiff_exogCO2), 4); % simulation x time x climsens_node x 1

% Next taking expectations over climate sensitivity ("cs" is
% for climate sensitivity)
margutility_cs = sum( bsxfun(@times,reshape(weights_climsens,[1 1 size(weights_climsens,1) 1]), margutility_dam), 3); % simulation x time x 1 x 1
margutility_cs_exogCO2 = sum( bsxfun(@times,reshape(weights_climsens,[1 1 size(weights_climsens,1) 1]), margutility_dam_exogCO2), 3); % simulation x time x 1 x 1
consdiff_cs = sum( bsxfun(@times,reshape(weights_climsens,[1 1 size(weights_climsens,1) 1]), consdiff_dam), 3); % simulation x time x 1 x 1
consdiff_cs_exogCO2 = sum( bsxfun(@times,reshape(weights_climsens,[1 1 size(weights_climsens,1) 1]), consdiff_dam_exogCO2), 3); % simulation x time x 1 x 1

% Obtain expected values at each time ("mc" is for Monte Carlo)
margutility_mc = mean(margutility_cs,1); % 1 x time x 1 x 1
margutility_mc_exogCO2 = mean(margutility_cs_exogCO2,1); % 1 x time x 1 x 1
consdiff_mc = mean(consdiff_cs,1); % 1 x time x 1 x 1
consdiff_mc_exogCO2 = mean(consdiff_cs_exogCO2,1); % 1 x time x 1 x 1

% Reshape to make dimensions be time x otherloop
margutility_mc = margutility_mc'; % time x 1
margutility_mc_exogCO2 = margutility_mc_exogCO2'; % time x 1
consdiff_mc = consdiff_mc'; % time x 1
consdiff_mc_exogCO2 = consdiff_mc_exogCO2'; % time x 1


% loop through times in which want scc
for scc_year = 1:Params.years_for_scc
       
    timevector = [scc_year:Params.horizon+(scc_year-1)]';
    
    % Convert utility to present value consumption, adjusting for population
    margutility_time = bsxfun(@times, margutility_mc(timevector,:), bsxfun(@times, exp(-Params.discount*(timevector-1)), (Params.cons_start*Paths.pop(1,timevector)'/Params.pop_start).^Params.power_cons ) );
    margutility_time_exogCO2 = bsxfun(@times, margutility_mc_exogCO2(timevector,:), bsxfun(@times, exp(-Params.discount*(timevector-1)), (Params.cons_start*Paths.pop(1,timevector)'/Params.pop_start).^Params.power_cons ) );
       
    % Integrate over time
    scc(scc_year,index_otherloop) = trapz(margutility_time,1);
    scc_exogCO2(scc_year,index_otherloop) = trapz(margutility_time_exogCO2,1);
    
end % end looping through scc_year

% convert units from Gt C to tCO2
scc(:,index_otherloop) = scc(:,index_otherloop)*1e-9/Params.co2_per_c;
scc_exogCO2(:,index_otherloop) = scc_exogCO2(:,index_otherloop)*1e-9/Params.co2_per_c;



%% Test effects of the perturbation

timetest = [1:100];
test_perturb(timetest,1,index_otherloop) = squeeze(max(max(Paths_perturb.M1(1,timetest,:,:) + Paths_perturb.M2(1,timetest,:,:) - Paths_perturb.M1_exogCO2(1,timetest,:,:) - Paths_perturb.M2_exogCO2(1,timetest,:,:),[],3),[],4));
test_perturb(timetest,2,index_otherloop) = squeeze(max(max(Paths_perturb.temp(1,timetest,:,:) - Paths_perturb.temp_exogCO2(1,timetest,:,:),[],3),[],4));
test_perturb(timetest,3,index_otherloop) = squeeze(max(max(Paths_perturb.cons(1,timetest,:,:) - Paths_perturb.cons_exogCO2(1,timetest,:,:),[],3),[],4));


%% Set up for next otherloop

if length(otherloop_vector)>1
    clear Paths Wiener Paths_perturb;
    clear nodes_climsens weights_climsens nodes_climsens nodes_damage weights_damage;
end

if index_otherloop == length(otherloop_vector)
    save(['asset_results.mat']);
else
    save(['asset_results' num2str(index_otherloop) '.mat']);
end
if index_otherloop > 1
    delete(['asset_results' num2str(index_otherloop-1) '.mat']);
end